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Abstract Orbit determination is possible for a chaotic orbit of a dynamical system, 
given a finite set of observations, provided the initial conditions are at the central time. 
In a simple discrete model, the standard map, we tackle the problem of chaotic orbit 
determination when observations extend beyond the predictability horizon. If the orbit 
is hyperbolic, a shadowing orbit is computed by the least squares orbit determination. 
We test both the convergence of the orbit determination iterative procedure and the 
behaviour of the uncertainties as a function of the maximum number n of map iterations 
observed. When the initial conditions belong to a chaotic orbit, the orbit determination 
is made impossible by numerical instability beyond a computability horizon, which 
can be approximately predicted by a simple formula. Moreover, the uncertainty of the 
results is sharply increased if a dynamical parameter is added to the initial conditions 
as parameter to be estimated. The uncertainty of the dynamical parameter decreases 
like n“ with a < 0 but not large (of the order of unity). If only the initial conditions 
are estimated, their uncertainty decreases exponentially with n. If they belong to a 
non-chaotic orbit the computational horizon is much larger, if it exists at all, and 
the decrease of the uncertainty is polynomial in all parameters, like n“ with a ~ 
1/2. The Shadowing Lemma does not dictate what the asymptotic behaviour of the 
uncertainties should be. These phenomena have significant implications, which remain 
to be studied, in practical problems of orbit determination involving chaos, such as 
the chaotic rotation state of a celestial body and a chaotic orbit of a planet-crossing 
asteroid undergoing many close approaches. 

Keywords Chaotic Motions • Numerical Methods 


1 Introduction 

Chaotic dynamical systems are characterized by the existence of a predictability horizon 
in time, beyond which the information on the state available from the initial conditions 
is not enough for meaningful predictions. Thus it appears a difficult task to perform 
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an orbit determination for a chaotic dynamical system, at least when the observations 
are spread over a time-span longer than the predictability horizon. 

Nevertheless there are practical problems of orbit determination in which the sys¬ 
tem is chaotic and the time-span of the observations is very long. It is important to 
understand the behaviour of the solutions, with their estimated uncertainties, in par¬ 
ticular when the variables to be solved for include not just the initial conditions but 
also some dynamical parameters. If the number of available observations grows, but 
simultaneously the time interval over which they are spread grows up to values com¬ 
parable to the predictability horizon, does the so lution become more accura te, and is 
the iterative procedure of differential corrections i Milani and Gronchil . 1201^ . Chap. 5) 
to find the least squares solution still possible? 

In this paper we use a model problem, namely the discrete dynamical system de¬ 
fined by the standard map of the pendulum, with just one dynamical parameter, the 
fj, coefficient appearing in equation ©• We also set up an observation process in which 
both coordinates of the standard map are observed after each iteration. In the obser¬ 
vations we include a simulated random noise with a normal distribution. Then, each 
experiment of orbit determination is also a concrete computation of a segment limited 
to n iterations (of the map and of its inverse) of a e shadowing orbit for the 5 pseudo 
trajectory defined by the observation process. The Shadowing Lemma (see Section El 
provides a mathematically rigorous result on the availability of shadowing orbits, but 
thanks to the orbit determination process we make explicit the relationship between e 
and 5 (see Section EH, and we explicitly compute the e-shadowing orbit. 

At the same time, each experiment provides an estimate of the standard deviation of 
each of the variables, including initial conditions and the parameter. These estimates 
can be plotted as a function of n, thus showing the relationship between accuracy, 
number of observations and time interval, measured in Lyapounov times (see SectionEJ. 

Of course the numerical experiments are limited to a finite number of iterations, 
while the Shadowing Lemma refers to an infinite orbit. However, the maximum number 
of iterations is controlled by another time limit, the computability horizon due to round 
off error. This limit can be estimated approximately by a simple formula, and it is found 
in numerical experiments as a function of both the initial conditions and the numeric 
precision used in the computations. 


1.1 Wisdom hypothesis 


In 1987 J. Wisdom was discussing the chaotic rotation state of Hyperion, when he 
claimed that numerical experiments indicated that the knowledge gained from measure¬ 
ments on a chaot ic dynamical sy stem grows exponentially with the timespan covered by 
the observations I WisdorJ . Il987l 'l. This pertained in particular to the information on 
dynamical parameters like the moments of inertia ratios for Hyperion, as well as the 
rotation state at the midpoint of the time interval covered by the observations, which 
he proposed would be determined with exponentially decreasing uncertainty. 

Therefore Wisdom suggests that the orbit determination for a chaotic system might 
be in fact more effective than for a non-chaotic one. It is clear from the context that he 
was referring to numerical results, thus his statement can only be verified with finite 
computations as close as possible to a realistic data processing of observations of a 
chaotic system with dynamical parameters to be determined. 












3 


We have set as a goal in this paper to test the behavior of the uncertainty in 
the dynamical parameter of our model problem. We shall discuss the implications for 
Wisdom’s claim in Section EU 


1.2 Application to planet-crossing asteroids 

In our solar system there are planet-crossing minor bodies, including asteroids and 
comets, by definition such that their orbits can, at some times, intersect the orbit of 
the major planets. In particular many of the Near Earth Asteroids (NBAs) can intersect 
the orbit of the Earth. These orbits are necessarily chaotic, at least over the timespan 
accessible to accurate numerical computations. 

Unfortunately, these orbits are especially important and necessary to be studied 
because of the very reason of chaos, namely close approaches to the major planets 
including the Earth: these approaches may, in some cases, be actual impacts on a finite 
size planet. 

The attempt to predict possibility of impacts by NEAs, in particular on our planet, 
is called Impact Monitoring, and it is indeed a form of orbit determination for chaotic 
orbits. There is a subset of cases of NEAs for which non-gravitational perturbations, 
such as the ones resulting from the Yarkovsky effect, are not negligible in terms of Im¬ 
pact Monitoring because of the exponential divergence of nearby orbits which amplifies 
i nitially very small perturbations _ 

( [Faxnocchia et 1 20131 : ICheslev et al.l . l2014l : lEarnocchia and Cheslevl . l2014l : ISooto et al.l . 

boidi . 

Thus the Impact Monitoring for these especially difficult cases is an instance of 
orbit determination of a chaotic system, with as parameters the 6 initial conditions 
and at least one dynamical parameter, such as a Yarkovsky effect coefficient to be 
solved for. We shall show in Section [5]2] that the weak determination of the dynamical 
parameter is a key feature of these cases. 


2 Orbit determination for the standard map 


The simplest example of a conservative dynamical system which has both chaotic and 
ordered orbits can be built by means of an area preserving map of a two dimensional 
manifold: 


Vk) 


/ ^k T Vk+l 

1 Vk+l = Vk -MsinXfe 


( 1 ) 


where fj. is the perturbation parameter, and S is the standard map. The system has 
more regular orbits for small fj,, and more chaotic orbits for large p.. We choose an 
intermediate value of p, e.g. p — 0.5, in such a way that ordered and chaotic orbits are 
both present. Figure [T] shows the strongly chaotic region around the hyperbolic fixed 
point, and a few regular orbits on both sides. 

The advantage of such example is that the least square parameter estimation pro¬ 
cess can be performed by means of an explicit formula. 

First we compute the linearized map 


DS = 


dxkj,.! dxk+1 
Xk Vk 

dyk+i dyk+i 
Xk Vk 


/ 1 - pcos{xf,) l\ 
\ -pcos{xk) ij 
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Standard map (mu=0.5) 



Fig. 1 Orbits of the standard map for the perturbation parameter /i = 0.5. Plotted is a blow 
up of the central region around the hyperbolic fixed point, showing the strongly chaotic region 
and a few regular orbits on both sides. 


and from this the linearized state transition matrix 

A ^ 9{Xk, Vk) 

^ d{xo, yo) 

which is the solution of the variational equation (for infinitesimal displacement in the 
initial conditions), and given by the recursion: 


^fe+l = DS Ak ; Aq = I . 

The variational equation for the derivatives with respect to the model parameter y is: 

d{xk+i, Vk+i) ^ djxk, Vk) dSf, 
dfi. dy dfi 

^ „ djxk, Vk) (- sin(a:fe) \ 

dy V-sin(a;fe)y 

Then we set up an observation process, in which both coordinates x and y are observed 
at each iteration, and the observations are Gaussian random variables with mean Xk 
(yk respectively) and standard deviation cr: we use the notation Xk{yo,o-) to indicate 
that the probability density function of the observation Xk is the normal A/’(xfe(/io), cr^) 
one, and similarly for yk- The residuals are: 


Cfc = Xk{yo,(x) - Xkiyi) 
ik = ykiyo,o-) - ykiyi)- 


( 2 ) 
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for k = —n ,..., n. In © ®/c(M0i <^) a-iid J/fc(M0i <^) a^re the observations at each iteration, 

fiQ is the “true” value and = /rg + is the current guess. _ 

_ Then the least squares fit is obtained from the normal equations ( Milani and Gronchil . 

bninl'i : 


C= ^ BjSfe ; D = -^ B_ 


k— — r 


T ( ik 
' ik 


Bu — 


9i^k,^k) 

d{xo,yo,fJ.) 


= - A 


k— — n ^ 

d{xk,yk) 

djj, 


( 3 ) 


The least squares solution for both, the parameter y and the initial conditions, is: 


r = c-^ 


with r the covariance matrix of the result. This is enough to find the least squares 
solution by iteration of the above differential eorrection. However, to assess the un¬ 
certainty of the result, weights should be assigned to the residuals consistently with 
the probabilistic model, in this case each residual needs to be divided by its standard 
deviation ct; then the distribution of the result {x,y,y) is a normal distribution with 
covariance matrix cr^ F. 


3 Shadowing Lemma 


The shadowing problem is that of finding a deterministic orbit as close as possible to a 
given noisy orbit. The so-called Shadowing Lemma is the main re sult about sh adowing 
near a hyperbolic set of a diffeomorphism. Anosov ( Anosov|.ll967ll and Bowen I BowerJ 
Il975l i proved the existence of arbitrarily long shadowing solutions for invertible hype r- 
bolic maps. Here we give an overview of these classical results, as in ( Pilvuginl . ll999l l. 

Let {X, d) be a metric space and let ^ be a homeomorphism mapping X onto itself. 
A 5-pseudotrajectory of the dynamical system ^ is a sequence of points Q = {x^. £ X : 
fe G Z} such that the following inequalities 


d{F{xk),Xk+i) <5. (4) 

hold. For a graphical description of a 5-pseudotrajectory, see Fig. [2] Usually, a 5- 



Fig. 2 A (5-pseudotrajectory. 


pseudotrajectory is considered as the result of using a numerical method to compute 
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Fig. 3 An £-shadowing. 


orbits of the dynamical system e.g., because of round off error. We say that a point 
X £ X shadows a pseudotrajectory C, = {xjt} if the inequalities 

d{<l>''{x),Xk) < e (5) 

hold (see Figure |3]|. If only one dynamical system ^ is considered, we will usually write 
e-shadows The existence of a shadowing point for a pseudotrajectory ( means that 
^ is close to a real trajectory of <P. 

The following statement is usually called the Shadowing Lemma. 

Theorem 1 If A is a hyperbolic set for a diffeomorphism !>, then there exists a neigh¬ 
borhood W of A such that for all e > 0 there exists 5 > 0 such that for any 5- 
pseudotrajectory with initial conditions f £ W there is a point x that e-shadows C,. 

The Anosov shadowing requires the existence of a hyperbolic set. It means that at 
each point there are two directions where the motion is either exponentially expanding 
(stable manifold) or exponentially contracting (unstable manifold). 

Definition 1 We say that a set A is hyperbolic for a diffeomorphism <I> G (^^(R") if: 

(a) A is compact and ^-invariant; 

(b) there exist constants C > 0, Aq G (0,1), and families of linear subspaces S{p), U(p) 
of R”, p G A, such that 

(b.l) S(p)©t/(p) = R"; 

(b.2) D${p)T{p) = T{${p)), pgA,T = S,U-, 

(b.3) 


|D^™(p)u| < C'Ao'lul for V G S{p), m > 0; 

\D<I>~^{p)v\ < CXo^\v\ for V G U{p), m > 0; 

The definition of a hyperbolic set is strictly related to the one of Lyapounov expo¬ 
nent: for each orbit inside a hyperbolic set, the Lyapounov exponents must be either 
> log(Ao) or < -log(Ao). 


3.1 Shadowing Lemma and orbit determination 

Our goal is to connect the Shadowing Lemma with the chaotic orbit determination, 
involving the least squares fit and the differential corrections. 

First of all we build a 5-pseudotrajectory by using a simulated observations pro¬ 
cess. In Section [2] we have set up such an observations process, with observations 
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{Xk{fj-o,o-),yk{fj-o,cr))- We claim that a sequence C = {{xk{yo,'x),yk{yo,(x))} is a S- 
pseudotrajectory for the dynamical system (a:oi 3/o)i with 5 = — yo\ + ICa, 

/C G R. To obtain this resnlt we compute the Euclidean distance; 

diSfj.* {xkiyo, o-),ykiyo,<x)), [xk+i (t^o, <x), Vk+i (yo, o-))) (6) 

For the sake of simplicity {xk+i,yk+l) ar® the observations, i.e. Gaussian random 
variables with mean Xk+i respectively), and standard deviation a, as in Sec. [21 

and Sfj,* {xk,yk) = (®fc+li i/fe+l)- Using these notations, (|6]) turns into 


d{S^j,*{xk,yk), ixk+i,yk+l)) = \J{Xk+l - Xk+iY + (i/fe+l - yk+lY 

We compute separately the two differences. 

liifc+i -yk+i \ = \yk+i - y* sinxk -yk+i -A/'(o,cr^)| 

= |A/’(0, 2(7^) — y* sin®/; cos(A/’(0, cr^)) + y* sin(A/’(0, cr^)) cos Xk + yo sini^l 
<N{Q,2a'^) + \yo-y\ ( 7 ) 

To allow the last estimate, we need to solve a technical problem: the Shadowing 
Lemma uses a uniform norm, that is the maximum of the distance between the 5- 
pseudotrajectory and the e-shadowing. On the contrary, the natural norm for the 
residuals of the fit is the Euclidean norm with the square root of the sum of the 
squares. However, since the number of residuals is not only finite but sharply limited 
by the numerical phenomena discussed in Sectional in a given experiment we can just 
take the maximum absolute value of the residuals and it is going to be /Cct, with 1C a 
number which in practice cannot be too large, e.g., in our experiment 1C = 5.9. 

Then we can approximate the quantities 0(a) and smaller, e.g., cos(A/'(0, cr^)) ~ 1 
and sin(A/'(0, CT^)) ~ 0. The x coordinate gives a similar result: 

l^fc+1 T yk+1 Xk 2/fc+l W*(0, (7 )| 

< A/’(0, cr^) +JC{0, 2(t^) + \yQ — y*\ = A/’(0, 3(t^) + \yo — y*\ (8) 

Putting together © and ((SJ we obtain 

\J (yfc+1 ?/fc+l)^ 

< V(A^(0, 3a2) + \yo - + (A7(0, 2a^) + \yo - y*\)^ < (9) 

< V2\yo — y*\ + \J (A/’(0, 3cr2))2 _|_ (W'(0, 2a’^)Y < '/2\yo — y*\ + JCa 

with /C G R. 

Therefore the sequence generated by the observations is a (5-pseudotrajectory for 
the dynamical system with S — y/2\yQ — /i*| -|- ICa. 

Figure [4] is an example of observations as a 5-pseudotrajectory. The observations 
are built with a = 10~^ and yo = 0.5, and the dynamical system is 5^*, with \y* — y\ = 
10”^; the circles have radius 5. 

The solution of the least squares fit (to the observations from —n to n), obtained 
by convergent differential corrections, is a finite e-shadowing, valid for the iterations 
from —n to n. 

We choose a value e > 0, that is a boundary on the maximum of the residuals. 
Then we choose a < ejK, and we set up the observations process. Next, we create 
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Delta-pseudotrajectory for the standard map 
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Fig. 4 An example of a (5-pseudotrajectory. Initial conditions are xq = 3, j/o = 0, fio = 0.5. 
S/i = 10“^, and a = 10“®. 


a first guess: a new orbit obtained with a small change of the initial conditions and 
of the dynamical parameter fj,: {{xk{ng),yi^{jj.g))} = Sfig{xg,yg), with Xg = xq + dx, 
Vg ~ yo~^dy, and y,g = fio + dji. Then we apply the differential corrections to the orbit. 
If the iterations converge, that is the last correction is very small, the maximum of the 
norm of the residuals is less than e (because the individual residuals are less than 3cr). 

At convergence, we obtain an initial condition {x*, y*) and a value of the dynamical 
parameter fi*, such that {x*,y*) is the (e, )-shadowing for the J-pseudotrajectory 
with 5 = V2\yt* — /tqI + ICa, for all the points used in the fit. 

The most important requirement is the convergence of the differential corrections, 
otherwise we cannot obtain the e-shadowing. This is far from trivial, because the chaotic 
divergence of the orbits introduces enormous nonlinear effects, for which the linearized 
approach of differential corrections may fail. To guarantee convergence, first we select 
the initial conditions xo,yo to be at the center of the observed interval, otherwise the 
initial conditions would be essentially undetermined along the stable manifold of the 
initial conditions. Second, we use a progressive solution approach, namely, given the 
solution with 2n-|-l observations with indexes between —n and n, we use the convergent 
solution XqjPq, fj,* for n as hrst guess for the solution with 2n-|-3 observations (between 
—n — 1 and n -I- 1). Then the initial guess is actually used only for the solution with 
3 observations, for which the nonlinearity is negligible. Still the convergence of the 
differential corrections depends critically upon the number n of iterations of the map, 
as explained in Section |4] 
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4 Numerical results 

The experiment was performed with the initial conditions at iq = 3 and i/o = 0, 
and the value of the dynamical parameter po = 0.5. The dynamical context for this 
orbit can be appreciated from Figure [T] showing that the initial conditions are indeed 
in a portion of the initial conditions space containing mostly chaotic orbits. For the 
observation noise we have used a standard deviation a — 


4.1 Computability horizon 


Determinant and eigenvalues of the state transition matrix 



Fig. 5 The eigenvalues and the determinant of the state transition matrix in a semilogarithmic 
scale, as a function of the number of iterations. Also shown is the linear fit to the large 
eigenvalue based on the first 180 iteration, with slope +0.091. The computation is in double 
precision and the number of iterations of the standard map n is 300 with the map and 300 
with its inverse. The determinant of the state transition matrix would be 1, for all n, in an 
exact computation. The numerical instability occurs when the eigenvalues reach the critical 
values sj'^lsd marked by the dotted lines. 


Figures and [6] show the absolute value of the eigenvalues of the state transition 
matrix forward and backward. The product of two eigenvalues should be 1 in exact 
arithmetic. When the condition number of the matrix becomes larger than the inverse 
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Determinant and eigenvalues of the state transition matrix 



Fig. 6 The eigenvalues and the determinant of the state transition matrix in a semilogarithmic 
scale, as a function of the number of iterations. Also shown is the linear fit to the large 
eigenvalue based on the first 300 iteration, with slope +0.086. The computation is in quadruple 
precision and the number of iterations of the standard map n is 800 with the map and 800 
with its inverse. The numerical instability occurs when the eigenvalues reach the critical values 
y/^, y/llsq marked by the dotted lines. 


of the machine rounding off error, the computation of the matrix becomes numerically 
impossible, and the computed value of the determinant is far from 1. 

In Figure[3the computations are performed in the standard double precision, that 
is with a mantissa of 52 binary digits and a round off relative error of s^i = = 

1.1 X 10”^®. We observe a numerical instability after 180 iterations: the determinant 
deviates from the exact value of 1 and the small eigenvalue starts increasing; the large 
eigenvalue keeps increasing, but there is a slight change of slope. Then we fit the slope 
of the large eigenvalue curve for the first 180 iterations, and get a Lyapounov indicator 
+0.091: it approximates the maximum Lyapounov exponent x the orbit to which 
our differential corrections converg^. 


^ There is no way to rigorously compute the Lyapounov exponents: in practice Lyapounov 
indicators extracted from finite propagations are used to assess, but not rigorously prove, the 
chaotic nature of the orbits. Note that it is a numerically well documented phenomenon that 
the indicators are not constant, but actually depend upon the time interval over which they 
are computed, although in most cases these changes are not very large and the conclusion that 
an orbit is chaotic is reliable. 
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The Lyapounov time is = 1/x, in this example ~ 11. To reach a ratio of 
eigenvalues of the state transition matrix of l/e^ we need a number of Lyapounov times 
log(l/y^), in this case ~ 18.4~ 202 iterations of the map. At about this number 
of iteratious the maximum and minimum eigenvalues of An are so widely apart in size 
that a bad conditioning horizon is reached, and the computation of the state transition 
matrix becomes numerically inaccurate. Hence near ±18.we observe the numerical 
instability in the computation of the determinant and of the eigenvalues of the state 
transition matrix. 

Figure [6] shows the same computations, with the same initial condition, but in 
quadruple precision, with a 112 bit mantissa and Eq = 2“^^^ = 9.6 x 10“^^. The 
change of slope in the eigenvalues curves occurs after ~ 300 iterations, while a full 
blown numerical instability occurs after ~ 550 iterations. The fit to the large eigenvalue 
for the first 300 iterations gives a Lyapounov indicator ±0.086, not very different from 
the one obtained in double precision. Thus we would expect the numerical instability 
to occur after log(l/y^)ri ~ 39.2 = 455 iterations. It appears that the rate of 

divergence decreases after 300 iterations, as shown by the change in slope, allowing to 
maintain at least the determinant near the exact value for about 100 more iterations 
beyond the value predicted above. 

The computability horizon represents the maximum number of iterations we can 
reach, before the computation becomes numerically unstable. The computability hori¬ 
zon strongly depends on the chaoticity of the system: more chaos, that is larger X) 
more instability; but also upon the precision of the computations. 

Thus, in the following we perform the numerical experiments in quadruple precision, 
in order to mitigate the problem of the numerical instability. We compute 500 iterations 
forwards and backwards, but we use only the first 300 iterations for the linear fits, to 
avoid the possibility that changes in slope, such as the ones apparent in Figure [6l 
contaminate our experimental results. 

The compatibly horizon is a hard limit in that it is not practically possible to 
increase the number of iterations to a much higher value. E.g., to push the horizon by 
a factor 10 above the value for double precision, we would need computations performed 
with real numbers represented with 800 byteO 

The conclusion is that the practical problem of chaotic orbit determination is mean¬ 
ingful only for a finite number of iterations, and the accuracy of the results can be tested 
only within the boundary of the computability horizon. 


4.2 Chaotic case 

Figure [3 shows the results in quadruple precision of the full 3-parameter fit: the 3 pa¬ 
rameters are the initial conditions and the dynamical parameter p. The determination 
of p is indeed not possible without simultaneous determination of the initial conditions. 

Even in quadruple precision we hnd a maximum value of n beyond which the itera¬ 
tive solution of the nonlinear least squares problem is divergent. This maximum turns 
out to be 599 in this experiment: it is close to what we have called the computability 
horizon, that is this limitation is due to the difficulty of computing the state transition 
matrix when the condition number is too large. 

^ Software to perform arithmetic computations with an arbitrary number of digits is avail¬ 
able, but the algorithms are too slow to be used even for our simple example. 
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Uncertainty 3 parameters 



Fig. 7 Standard deviation of the solutions for the initial conditions and for the dynamical 
parameter fi (continuous lines), and actual error (nominal solution of the fit minus real value 
used in the simulation, dashed lines), as a function of the number of iterations. 


The curves in Fig. [7] represent both the formal standard deviation and the actual 
error of the solutions of the least squares fit, as a function of n in a semilogarithmic 
plot. Both the formal standard deviation and the actual error of fj, do not decrease 
exponentially. Indeed, in Fig. [S] we have the same behavior of the curves that we have 
already seen in Fig. 0 but in a log-log plot, in which a constant slope a would imply 
a power law proportional to n“. The slopes of the lines that fit the uncertainties are: 
—0.675 for the dynamical parameter jj,, —0.833 and —12.030 for the initial conditions x 
and y, respectively. This plot in logarithmic scale shows that the uncertainty for y and 
X does not decrease exponentially. It is also apparent that one of the initial conditions 
(y) is better determined than the other one (x), with an improvement as a function 
of n which could be exponential. This is a property of the specific initial condition we 
have used, for other choices we can get three parameters determined with comparable 
accuracy, none of them with exponential improvemen10. 

Figure [9] shows the results for the standard deviation and the actual error when 
solving only for the initial condition. The 2x2 portion of the normal matrix which refers 
only to the initial conditions is not badly conditioned. Also as a result of this, we are 
able to get convergence of the differential corrections up to ±742 iterations, which is 

® This depends upon the orientation of the stable and unstable directions at the initial 
condition. 
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Decrease of uncertainty and fit (3 par: x, y, p) 
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Fig. 8 Uncertainty of the solution of the least squares fit for the initial conditions and for the 
dynamical parameter /r in a logarithmic scale. 

even beyond the numerical stability boundary. If the fit is done by using only up to 
300 iterates, to avoid the apparent slope change, the slopes shown in this Figure are 
—0.084 for X and —0.083 for y, note that the Lyapounov indicator for the same interval 
is +0.086. Thus exponentially improving determination of the initial conditions only is 
possible, and the exponent appears to be very close to the opposite of the Lyapounov 
exponent. 


4.3 Ordered case 

An ordered case can be obtained with a change of the initial conditions. For the nu¬ 
merical experiments we have chosen xg = 2, yo — 0 and pg = 0.5. In the ordered case 
we have not the problem of the computability horizon and the Lyapounov exponent 
is very small: actually, it could be zero if we are on a Moser invariant curve. Thus we 
have computed 5000 iterations. 

Figure [TU] gives a summary of our numerical experiment in the ordered case. The 
Lyapounov indicator is very small (~ 10~^), and can be made even smaller by continu¬ 
ing the experiment for larger values of n. As a consequence, the state transition matrix 
is not badly conditioned, and the computability horizon is much beyond the number 
of iterations we have used (if it exists at all). Thus the lack of chaoticity implies the 
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Uncertainty 2 parameters 



Fig. 9 Standard deviation of the solutions for the initial conditions (continuous lines), true 
errors for the same 2 parameters (dashed lines). 


practical absence of the computability horizon, and we can determine all the param¬ 
eters with very good accuracy, even if we are not in exact arithmetic. The values of 
the slopes of the fit to the uncertainty are —0.504 for —0.504 and —0.488 for x and 
y respectively, the corresponding regression lines are shown in the log-log plot on the 
bottom right. As it is clear by comparing the top right and the bottom left plot, the 
standard deviation for the solution with only 2 parameters have very much the same 
behavior, indeed in a log-log plot (not shown) we can get slopes —0.511 for x, —0.481 
for y. 

All these power laws are close to the inverse square root of the number of iterations, 
namely the same rule as the standard deviation in the computation of a mean. We do 
not have a formal proof of this, but we conjecture that for an orbit on a Moser invariant 
curve (for which the Lyapounov exponents are exactly zero) the standard deviations 
for all the parameters decrease as 1/ y/n. 


5 Conclusions 

We have understood the concept of computability horizon as a consequence of numer¬ 
ical instability in the computation of the state transition matrices, thus providing a 
comparatively simple empirical formula to approximately predict the horizon. This is 
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Determinant and eigenvalues of the state transition matrix 



Uncertainty 2 parameters 




Decrease uncertainty and fit (3 par: x, y, p) 



Fig. 10 Top left: eigenvalues of the state transition matrices, for the chosen regular initial 
conditions and for ±5000 iterations. Top right: solutions for the initial condition only, from 
the top condition number of the normal matrix, standard deviation of y and for x. Bottom 
left: solutions for three parameters, from the top condition number, standard deviation of x, 
of y, of /i. Bottom right: log-log plot of the 3 standard deviations, with very similar slopes. 


a practical limitation which applies to any attempt at orbit determination of chaotic 
orbits. 

We have used numerical experiments in quadruple precision, but nevertheless lim¬ 
ited by the computability horizon to few hundreds of iterations for the very chaotic 
orbit used in our test. From these we have found the following three empirical facts: 

1. If only initial conditions are determined for a chaotic orbit, the uncertainty can 
decrease exponentially with the number n of iterations of the map, and the exponent 
of this decrease is close to a Lyapounov exponent. 

2. If a dynamical parameter jj. is determined together with the initial conditions of a 
chaotic orbit, the decrease in uncertainty is polynomial in n for p and for at least 
one of the initial coordinates. 

3. If the initial conditions belong to an apparently ordered orbit, that is such that 
there is no evidence of a positive Lyapounov exponent, it is possible to determine si¬ 
multaneously /i and the initial conditions with uncertainty decreasing polynomially 
with n. The case in which only the initial conditions are determined gives the same 
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result. Moreover, all the power laws n“ for these uncertainties appear numerically 
to have a ~ —1/2. 


5.1 Discussion on the Wisdom hypothesis 

The statement by Wisdom, as a practical rule for concrete orbit determination, ap¬ 
pears to be first limited by the computability horizon. Second, the actual decrease of 
the uncertainty, going as far as it can be done numerically, is not exponential, but 
polynomial, as n“, with a negative and rather small, although we have found that the 
value of a depends upon the initial condition^. Note that the orbit determinations in 
which the only parameters to be solved are the 2 initial coordinates show an exponen¬ 
tial decrease as exp(—an), where a appears to be close to the Lyapounov exponent x, 
but the strong correlations appearing when 3 parameters are solved degrade the result 
in a very substantial way. 

This needs to be compared to the regular case, shown in Figure 1101 where the 
standard deviations for each of the 3 fit parameters decrease approximately according 
to an 1/%/n law, as prescribed by the standard rule for the estimate of the mean with 
errors having a normal distribution. Indeed it is possible that the determination of fi for 
some chaotic cases, including the example shown in Figure [8l decreases faster than for 
an ordered case, but the decrease is anyway polynomial, proportional to n“ with some 
different negative a, thus the difference is not very large, given the tight constraint on 
the maximum possible value of n. 


5.2 Examples from Impact Monitoring 


One feature of our results is that adding a dynamical parameter to the list of parameters 
to be determined results in degradation in the normal matrix, thus in much slower 
decrease of the uncertainties as the number of observations grows. The problems of 
orbit determination for NEA undergoing several close approaches to the Earth (or 
other planets) is more complex than our simple model, but we have found that the 
phenomenon described above does occur in a remarkably similar way. 

In Figure fm we show two probability distributions, as derived from the orbit de¬ 
termination of the asteroid (410777) 2009 FD. The narrow peaked distribution corre¬ 
sponds to an orbit determination with 6 parameters, the initial conditions only: the 
standard deviation is 6 x 10^ km. The much wider distribution corresponds to a fit 
with 7 parameters, including the constant A 2 appearing in the transverse acceleration 
due to the Yarkovsky effect: the STD is ~ 2.3 x 10® km. The Yarkovsky effect is a form 
of non-gravitational perturbation due to thermal radiation emitted anisotropically by 
the asteroid, and is indeed very small. However, when the uncertainty resulting from 
the covariance matrix of the orbit determination is propagated for ~ 170 years after 
the last observation available, not only the Yarkovsky effect has a long enough time to 
accumulate but it is also enhanced by th e exponentia l diver gence of nearby orbits, the 
Lyapounov time being about 15.3 years ( Snoto et al.l . Eoiih [Figure 5]. 


^ We are showing figures and giving data only for one initial condition, but of course we 
have run many tests. 
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Fig. 11 Two different Probability Density Functions (PDF) for the trace of possible solutions 
on the Target Plane of the close approach of asteroid (410777) 2009 FD to the Earth in the 
year 2185. Superimposed and on a different vertical scale are the keyholes relative to impacts 
in different years between 2185 and 2196; the height of the bar is proportional to the width 
of the keyhole, thus the Impact Probability can be computed as product of the width and the 
PDF. 


The practical consequence of this increase of the uncertainty arises from the fact 
that the Target Plane of 2009 FD for 2815 includes some keyholes, small portions cor¬ 
responding to impacts with the Earth (either at that time or a few years later, until 
2196). With the 7 parameters solutions these keyholes are within the range of outcomes 
with a significant value of the Probability Density Function, thus the impacts have a 
non-negligible probability, the largest being an Impact Probability of ~ 1/370 for 2185. 
If on the contrary the orbit was estimated with 6 parameters only, then the probability 
would appear to be even larger for an impact in 2190, and all the other keyholes (in¬ 
cluding the one for 2185) would correspond to negligible impact probabilities. Given 
that the impact, if it was to occur, would release an energy equivalent to 3, 700 Mega- 
Tons of TNT, this difference is practically relevant. In fact, the solution including the 
Yarkovsky effect leads to a more reliable estimate of the Impact Probabilities, because 
the Yarkovsky effect exists and needs to be taken into account. 

Is the discrepancy in the uncertainties with and without the dynamical parameter 
in the fit essentially the same phenomenon we have found in our simple model? We do 
not know the answer to this question, but we shall investigate this issue in the future. 
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